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Event simulations in a transport model for intermediate energy heavy ion collisions: 

Applications to multiplicity distributions 
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We perform transport model calculations for central collisions of mass 120 on mass 120 at labo¬ 
ratory beam energy in the range 20 MeV/nucleon to 200 MeV/nucleon. A simplified yet accurate 
method allows calculation of fluctuations in systems much larger than what was considered feasible 
in a well-known and already existing model. The calculations produce clusters. The distribution of 
clusters is remarkably similar to that obtained in equilibrium statistical model. 
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I. INTRODUCTION 


II. THE PRESCRIPTION 


There is an enormous amount of experimental and the¬ 
oretical work on multifragmentation in heavy ion colli¬ 
sions at intermediate energy. There are two classes of the¬ 
oretical models: (a) dynamical and (b) statistical where 
one assumes that because of two-body collisions the col¬ 
liding systems equilibrate and break up into many frag¬ 
ments according to availability of phase space. Unques¬ 
tionably the most used and popular dynamical model 
is the Boltzmann-Uehling-Uhlenbeck (BUU) model (also 
called by various other names). Here we restrict ourselves 
to BUU [J] . In the original formulation, the BUU model 
gave an account of one body properties @ and thus was 
not suitable to describe multifragmantation. Later it was 
extended to include fluctuations which made it suitable 
for event by event simulation 0. This is the focus of 
our attention. In the past event by event simulation in 
this model was limited to about mass number 30 on 30. 
The problem was a practical one, namely, it required a 
large computing effort. We show that with a slight re¬ 
formulation without changing any physics or numerical 
accuracy we can very significantly reduce the execution 
time and we can handle much larger systems. Computa¬ 
tion becomes as short as an ordinary BUU calculation. 
It is instructive to do large systems (finite number effects 
often hide important bulk effects) and more importantly, 
the fragmentation must be investigated over an energy 
range to unravel many interesting effects. The objective 
of doing examples is to demonstrate that many revealing 
features are seen. These also allow us to relate with other 
models. 

There are many models for multifragmentation. There 
are some which can be labeled as “quantum molecular 
dynamics” type 0,1 These are different in spirit to the 
model used here. Closer in spirit yet quite distinct are 
some studies based on a Langevin model mi where 
we have mentioned only a few. We will have occasion to 
refer very briefly to only a small number of these. The 
literature in the Langevin approach is huge. 


The basic features of our transport model calcula¬ 
tion are contained in the Boltzmann-Uehling-Uhlenbeck 
(BUU) model as developed in [lj and Q but some mod¬ 
ifications were made. For brevity we will almost entirely 
skip the physical motivations and details for the mod¬ 
els of refs, [lj and [3j as they are not only adequately 
discussed in the original papers but also elsewhere 0, Q 
where some different models are also introduced. The 
modifications we make to refs. 0 and 0 are discussed 
fully. 

The start of our consideration is the cascade model 
n Here each nucleus is considered as a collection of 
point nucleons whose positions are assigned by Monte- 
Carlo sampling. The projectile nucleus A approaches the 
target B with a beam velocity and two body collisions be¬ 
tween the nucleons take place. When these finish we have 
one event. We only consider B same as A and central col¬ 
lisions. It is convenient to run several events simultane¬ 
ously. Let us label the number of runs by N. In cascade 
the different runs do not communicate with each other. 
Thus nucleus 1 hits nucleus 1’, nucleus 2 hits nucleus 
2’, .•••nucleus N hits nucleus N\ In BUU we introduce 
communication between runs. What we were calling nu¬ 
cleons we now call test particles (abbreviated from now 
on as tp). The density p{f) is given by n/(Sr) 3 N where 
n is the number of tps in a small volume (Sr) 3 . As far 
as collisions go, in usual applications of BUU one still 
segregates different runs. By segregating the collisions 
one is able to use a nn , the nucleon-nucleon cross-section 
and reduce computation. If we considered collisions be¬ 
tween all tp’ s, the collision cross-sections would have to 
be reduced. In between collisions tps move in a mean 
field (Vlasov propagation). Applications of BUU as sum¬ 
marised above have met much success in explaining av¬ 
erage properties such as average collective flows etc. 

To explain multifragmentation, multiplicity n a as a 
function of a where a is the mass number of the com¬ 
posite, one needs an event by event computation in the 
transport model. Bauer et al made the following pre¬ 
scription 0. Now all tps are allowed to collide with one 
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another with a cross-section of a nn /N. Collisions are fur¬ 
ther suppressed by a factor N but when two tps collide 
not only those two but 2 (IV — 1) tps contiguous in phase 
space change momenta also. Physically it represents two 
actual particles colliding. When collisions cease we have 
one event. A second event needs a new Monte-Carlo of 
tps and then the evolution in time. 

The prescription we use here is the following. This is 
the middle ground between ref. [l] and ref. Q. As in 
ref. [J] for nucleon-nucleon collisions we consider 1 on 
l’(eventl), 2 on 2’(event2) etc. with cross-section a nn . 
For event 1 we will consider nn collisions only between 1 
and 1’. The collision is checked for Pauli blocking as in 
ref. HJ. If a collision between i and j in event 1 is al¬ 
lowed we follow ref. Q and pick N -1 tps from all the tps 
closest to i and give them the same momentum change 
A p as ascribed to i. Similarly we pick N-l tp’s closest 
to j and ascribe them the momentum change —Ap, the 
same as suffered by j. We will return to more details 
about this later. As a function of time this is continued 
till event 1 is over. For event 2 we return to time t= 0, 
the original situation (or a new Monte-carlo sampling for 
the original nuclei), follow the above procedure but con¬ 
sider nn collisions only between 2 and 2’. This can be 
repeated for as many events as one needs to build up 
enough statistics.The advantage of this over that used in 
ref. [3| is that here, for one event, nn collisions need to 
be considered between (Na+Nb) nucleons (A/ 4 =number 
of nucleons in A , IVb =number of nucleons in B ) whereas 
in method of ref. Q, collisions need to be checked be¬ 
tween (Na + Nb) x N tps. Hence, in our calculation, 
total number of combinations for two-body collision is 
reduced by a factor of 1/N 2 . Since typically N is of the 
order of 100 this is a huge saving in computation and has 
allowed us to treat mass as large as 120 on 120 over a 
substantial energy range. It is expected that the model 
used in ref. Q and the one used here will give similar 
results. The number of collisions for one event should be 
about the same in both prescriptions. The characteris¬ 
tics of scattering are the same. The objects that collide 
in our calculation arise from a coarse grain representa¬ 
tion of the initial phase space population of two nuclei. 
In ref. ff| a fine grain representation is used. But since 
many events are generated any difference should disap¬ 
pear. The Vlasov propagation is the same. For mass 
40 on 40 we compare our results with those using the 
method of ref. [3| (Fig.l). The agreement between the 
two calculations for multiplicities is remarkable. We re¬ 
gard our method as a very convenient short cut to the 
numerical modeling of ref. Q. The theoretical formula¬ 
tion in ref. [3j is more appealing and “democratic” but 
numerically our method gives indistinguishable results. 

One bonus of our prescription is that one sees some 
common ground between the BUU approach and the 
“quantum molecular dynamics” approach. In the latter 
nucleons are represented by Gaussians in phase space; 
the centroids have an r and a p which are originally gen¬ 
erated by Monte-Carlo calculations. These collide. This 



FIG. 1: (Color online) Comparison of mass distribution cal¬ 
culated according to the prescription of ref. Q (blue dotted 
lines) and the present work (red solid lines). The average 
value of 5 mass units is shown. The cases are for central col¬ 
lision of mass 40 on mass 40 for different beam energies (a) 
25, (b) 50 and (c) 100 MeV/nucleon. 500 events were chosen 
at each energy. 


corresponds to “nucleons” colliding in our prescription. 
As the centroids move after collision, they drag the Gaus¬ 
sians along. The Gaussian wave packets in position and 
momentum space provide the mean-field and Pauli block¬ 
ing. The Gaussians do not change their shapes or widths. 
These are very strong restriction and lead to very differ¬ 
ent mean field propagation. The Vlasov propagation has 
much more flexibility and originates from more funda¬ 
mental theory. 


III. SOME DETAILS OF THE SIMULATIONS 


We provide some details of the calculation. For Vlasov 
propagation we use the lattice Hamiltonian method jllj 
of Lenk and Pandharipande which accurately conserves 
energy and momentum. The mean field Hamiltonian for 
Vlasov propagation is also adopted from that work. The 
potential energy density is: 

p(r) 


The values of the constants are A=-2230.0 MeV/m 3 , 
B= 2577.85 MeV/m 7 / 2 , a= 7/6, /9 O =0.16/to" 3 , c=-6.5 
MeV/m 5 / 2 . The last term in the right hand side of eq.(l) 
gives rise to surface energy in finite nuclei. That favours 
the formation of larger composites, for example, the oc¬ 
currence of a nucleus of A nucleons over the formation 
of two nuclei of A/2 nucleons. Entropy works the other 
way. 

Further details are : 

(1) Calculations were done in a 200 x 200 x 200 /to 3 box. 
The configuration space was divided into 1 /to 3 boxes. 

(2) For results shown here the code was run from t = 
0/to/c to t = 200 fm/c. Positions and momenta of tps 
were updated every At = 0.3 fm/c. 
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(3) For nucleon-nucleon collision we follow Appendix. 
B. of ref. pj. 

(4) The number N was set at 100. 

(5) Once the two-body collisions are nearly over, con¬ 
tiguous boxes with tps that propagate together for a long 
time are considered to be part of the same cluster. The 
contiguous boxes have at least one common surface and 
the nuclear density exceeds a minimum value ( dmin ). Dif¬ 
ferent d m i n values as 0.002, 0.005, 0.01, 0.015 and 0.02 
fm -3 are tried to check the sensitivity of this parameter. 
It is observed that the fragment multiplicity distribution 
is not changing very much with d m i n , therefore we use 
dmin =0.01 fm -3 for further calculations. 


IV. RESULTS 

In Fig.2 we show plots of multiplicity against mass 
number a for 120 on 120. Four beam energies are shown. 
For each energy 1000 events were taken. We show results 
of averages for groups of five consecutive mass numbers. 
The prominent feature we wish to point out is that at 
low beam energy (50 MeV/nucleon) the multiplicity first 
falls with mass number a, reaches a minimum, then rises, 
reaches a maximum before disappearing. As the beam 
energy increases the height of the second maximum de¬ 



FIG. 2: (Color online) Mass distribution from BUU model cal¬ 
culation for Na = 120 on Nb = 120 reaction at beam energies 
(a)50 MeV/nucleon, (b)75 MeV/nucleon (c)100 MeV/nucleon 
and (d)150 MeV/nucleon. The average value of 5 mass 
units are shown. At each energy 1000 events are chosen. 
Only central collisions are considered here but even at E v = 
50 MeV/nucleon, nucleons in the peripheral region passes 
through and largest fragment remaining is less than the sum 
of the masses of the two nuclei. 



FIG. 3: Mass distribution from the Canonical Thermodynam¬ 
ical Model (CTM) calculation for fragmentation of a system 
of mass Ao = 192 at temperature (a)6.5 MeV, (b)7.5 MeV 
(c)10 MeV and (d)14 MeV. 


creases. At beam energy 75 MeV/nucleon the second 
maximum is still there but barely. At higher energy 
the multiplicity is monotonically decreasing, the slope 
becoming steeper as the beam energy increases. This 
evolution of shape of the multiplicity distribution is of 
significance as we will emphasize soon but let us point 
out this evolution of shape was a long time prediction of 
the canonical thermodynamic model (CTM) [HI,EH- For 
transport models the natural variable is the beam energy. 
For CTM the natural variable is the temperature T. For 
illustration we have shown the multiplicity distribution 
for a system of 192 particles in CTM at temperatures of 
6.5 MeV, 7.5 MeV, 10 MeV and 14 MeV (Fig. 3). The 
calculations with BUU and CTM are so different that 
the similarity in the evolution of the shape in multiplicity 
distribution is very striking. Indeed this correspondence 
provides the support for assumptions of statistical model 
from a microscopic calculation. 

To proceed further with the correspondence between 
the two models we need to establish a connection be¬ 
tween E p of BUU and temperature T of CTM. Temper¬ 
ature T of CTM will give an average excitation energy 
E* of the multifragmenting system in its center of mass 
(cm) 0 - We can calculate the excitation energy (E*) 
in the cm from (E p ) by direct kinematics by assuming 
that the projectile and the target fuse together. In that 
case the excitation energy is E* = A p E p /(A p +A t ) where 
A p and A t are projectile and target masses respectively. 
This value is too high as a measure of the excitation en¬ 
ergy of the system which multifragments. The nucleons 
at the edges of the two nuclei pass through carrying a lot 
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FIG. 4: (Color online) Top left curve (a) is a Canonical Thermodynamical Model (CTM) calculation for excitation (E*) vs. 
temperature(T) for Ao = 192. Between 6 MeV and 7.5 MeV temperatures, E* rises quickly. The dE*/dT slope increases 
sharply with mass size Ao and is indicative of first order phase transition. Bottom left curve (b) is also a CTM curve showing 
that the size of largest cluster drops sharply between 6 MeV and 7.5 MeV. Again this is first order liquid gas phase transition. 
Top right (c) is also with CTM but Amax/Ao is plotted against excitation energy per nucleon instead of temperature. The 
change of liquid to gas is necessarily slower, the range of energy for the change is dictated by latent heat. Bottom right (d) is 
the calculation from transport model. 


of energy and are not part of the multifragmenting sys¬ 
tem. These are pre-equilibrium particles. In experiments 
about 20% of nucleons are emitted as fast pre-equilibrium 
particles: see for example mm- Further details can be 
found in ref. but this is what we do basically. We 
go to the cm of the two ions to do BUU and at the end 
discard 20% particles (these have the highest energies), 
and measure energy (potential plus kinetic) of the rest. 
To find the excitation energy we subtract Thomas-Fermi 
ground state energy of the rest with the Hamiltonian of 
eq.(l). 

Fig.4 gives some CTM results and also makes a com¬ 
parison of one CTM result with transport model result. 
The top left diagram is E* vs T in CTM for 192 parti¬ 
cles (A 0 =192=80% of 240). This approximates usual E* 
vs T for first order phase transition. There is a boiling 
point temperature T which remains constant as energy 
increases. In the example here because we have a very 
finite system, the slope dE*/dT is not infinite but high. 
Let us now consider lower left diagram again drawn in 
CTM. Here A max is the average value of the largest clus¬ 


ter. A high value of A max /Ao means liquid phase and low 
values means gas phase. The criteria of deciding which 
composites belong to the gas phase and which to the liq¬ 
uid phase are discussed in detail in two previous papers 

gi (a. 

In the bottom left diagram, one sees more dramati¬ 
cally that in a short temperature interval liquid has trans¬ 
formed into gas. The only input in our transport model 
is the beam energy. The common dynamical variable in 
both our model and CTM is E*. Of course the E* in 
CTM is an average whereas the E* in transport model 
is a microcanonical E*. In the top right corner of Fig.4 
is the plot of A max /Ao as a function of E* in CTM. 
The transformation from liquid to gas is more gradual, 
essentially spanning the energy range across which, liq¬ 
uid transforms totally into gas. Even for a large system, 
where the transformation of liquid to gas as a function 
of temperature is very abrupt, the transformation as a 
function of energy per particle will be quite smooth. The 
bottom right in Fig.4 is from our transport model calcula¬ 
tion. The similarity with the CTM graph is close enough 







5 


that we conclude the transport model calculation gives 
evidence of liquid-gas phase transition. To find closer cor¬ 
respondence between transport model calculations and 
CTM, it will be best if we can deduce at least an approx¬ 
imate value of temperature for each beam energy. For an 
interacting system this is very non-trivial. Formulae like 
= ^f are obviously inappropriate. One might try to 
exploit the thermodynamic identity T = (^rf)v- This re¬ 
quires obtaining a value of the entropy for an interacting 
system. We will be working on this in future. 

In concluding this section, we mention that while 
we have established a correspondence between transport 
model results and CTM results, a more natural choice 
would have been to compare transport model results 
with multiplicity distributions obtained from the mi- 
crocanonical statistical multifragmentation model(SMM) 
[2Cj. These are not available to us. However for the only 
cases investigated we found that CTM and SMM results 
were quite close jU| so the correspondence we have found 
here between transport model results and CTM will pre¬ 
sumably hold for SMM also. Multiplicity distributions in 
16 O+ 80 Br were done with dynamic models before. The 
work in Ref. [22| used molecular dynamics. The work 
in Ref. [23] comes closer in spirit to ours. The colliding 
system was small and no attempt was made to link the 
work with statistical models or phase transition. 


V. DISCUSSION 

We now look at one feature of the model that raised 
concerns and led to a lot of work to propose alternative 
methods for calculations 0E3. This is related to dan¬ 
gers of crossing fermionic occupation limits in the model 
here (as in the model of ref. Q.). As mentioned already, 
if Pauli blocking allows two tps i and j to collide, then 
not only these two but also N — ltps closest to i and N -1 
closest to j move to represent that two actual nucleons 
scatter. The tps that move with i are denoted by i s . with 
s=0 to N — l.The square of the distance is taken to be 
dn s = (r *°p‘i is) + < ' Pio ~,[ is - 1 . Here R is the radius of the 

-tx p p 

static nucleus of A=120 and pf the Fermi momentum.. 
The tps j s are then chosen from the rest of the tps. Define 
now < pi >= , similarly < pj >. One then consid¬ 

ers a collision between < pi > and < pj > and obtain a 
A p for < pi > and - A p for < pj >. This A p is added to 
all p’is and — A p to all pj. Since the tps are moved with¬ 
out verifying Pauli blocking there may be cases where 
one exceeds the occupation limits for fermions. 

Initially the two ions have a very compact occupation 
at two different corners of phase space. Collisions make 
a far wider region of phase space available to nucleons so 
this problem may not be severe. An accurate estimation 
of exceeding the fermionic limit of occupation at vari¬ 
ous parts of phase space is very hard to compute in our 
present problem but some measures are relevant. 

For 120 on 120 at 100 MeV/nucleon beam energy (50 



Time (fm/c) 

FIG. 5: (Color online) Variation of average availability fac¬ 
tor (see text) with time (red line) for Na = 120 on Nb = 120 
reaction at beam energy 100 MeV/nucleon. The lower curve 
(black dotted line)is the average availability factor (/) at the 
phase space points of arbitrarily chosen 120 test particles in 
an isolated mass 120 nucleus as they move in time. The fluc¬ 
tuations from the value 0 reflects uncertainties, probably due 
to fluctuation in initial Monte-Carlo simulations. 


MeV/nucleon beam energy was studied also) we follow 
one event as a function of time. In every collision in the 
event, 2 N tps change momenta. To be specific, let the tp 
i s move from f,pT n to r,p~fi. We check if by moving i s to 
the final phase space point ( r,pji ) we cross the fermionic 
limit. We build a six dimensional unit cell in phase space 
around this final phase space point UM- The volume of 
the unit cell Should be small so that one is investigating 
the phase-space occupation very near (r,p/j), but it can 
not be too small since Monte-Carlo simulation has noise 
which can cloud the actual effect. In accordance with 
past calculations we chose a volume of unit cell in phase 
space where 8 tps is the maximum number allowed for 
fermions. If n is the number of tps (not including i s ) 
already in the unit cell we define an availability factor / = 
1.0 — n/8. If /=0 we are already at the limit of fermionic 
occupation. If / is negative we have crossed the quantum 
limit and are in the classical regime. Any positive number 
between 0 and 1 will accommodate additional fermion. 
For each collision there are 200 /’s to be calculated so 
for each collision we get an ave / and that is plotted in 
Fig.5. We have shown results for t=25 to 125 fm/c when 
most of the action takes place. For reference we also plot 
average / for randomly chosen 120 tps in a static mass 
120 nucleus as they move around in time. This number 
should ideally be 0 and not fluctuate. The deviations 
from zero in the static case probably largely arise due 
to fluctuations in Monte-Carlo sampling. This degree of 
uncertainty must be also present in the values of / we 
have plotted for collisions. In spite of these uncertainties 
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the predominantly positive values of / as displayed in 
Fig.5 lead us to believe that the general trends we find 
in our calculation will hold. 

If in a collision all of the 200 tp moved to locations 
where / were all positive we will stay within the fermionic 
limits. In case there is a tp which does not satisfy this 
we can try to improve the situation by discarding that 
tp and choosing the next available tp to be part of the 
cloud. Complications arise because when some of the 
previously chosen tps are discarded for new ones the av¬ 
erage momentum of the clouds will change, new A p will 
have to be used so the final resting spots obeying energy 
and momentum conservation will change too. An itera¬ 
tive procedure needs to be formulated but convergence 
may be slow. 

Alternative methods have been proposed. The two pa¬ 
pers which give procedural details of moving two clouds 
of tps from initial positions to final positions with a 
stricter adherence to fermionic limits are refs. 00 . 
Multiplicity distributions are not given so we can not 
compare. Even if the multiplicity distributions turn out 
to be similar, higher order correlations can be very differ¬ 
ent. The present work extended the first proposed model 
of fluctuations in BUU to a larger system at many ener¬ 
gies and a very interesting lesson was learned. The gross 
features of multiplicity distribution do resemble strongly 
the results from equilibrium statistical models which have 
proven very successful in explaining experimental data. 


VI. SUMMARY 

An event by event simulation of a transport model was 
made at collisions of moderately heavy ions at zero im¬ 
pact parameter. Multiplicity distributions were calcu¬ 
lated. They are remarkably similar to those obtained 
from a equilibrium statistical model (CTM). This work 
therefore justifies the use of the equilibrium statistical 
model for data fitting. This statistical model implies first 
order phase transition in large nuclear systems at finite 
temperature. It will be of interest to quantify more pre¬ 
cisely the correspondence of transport model result and 
statistical model results. That work is in progress. 
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